subroutine initcom

!----------------------------------------------------------------------- 
! 
! Purpose: 
! Initialize Model commons, including COMCON, COMHYB, COMMAP, COMSPE,
! and COMTRCNM
! 
! Method: 
! 
! Author: 
! Original version:  CCM1
! Standardized:      L. Bath, Jun 1992
!                    L. Buja, Feb 1996
!
!-----------------------------------------------------------------------
!
! $Id$
! $Author$
!
!-----------------------------------------------------------------------
   use shr_kind_mod,    only: r8 => shr_kind_r8
   use pmgrid,          only: plat, plev, plon, plevp
   use spmd_utils,      only: masterproc, iam
   use pspect
   use comspe
   use rgrid,           only: nlon, beglatpair, wnummax, nmmax, fullgrid
   use scanslt,         only: nlonex, platd, j1
   use gauaw_mod,       only: gauaw
   use commap
   use physconst,       only: rair, rearth, ra
   use time_manager,    only: get_step_size
   use abortutils,      only: endrun
   use scamMod,         only: scmlat,scmlon,single_column
   use cam_control_mod, only: nsrest
   use hycoef,          only: hycoef_init, hypd, hypm
   use eul_control_mod, only: ifax, trig,eul_nsplit
   use cam_logfile,     only: iulog
!-----------------------------------------------------------------------
   implicit none
!-----------------------------------------------------------------------

! Local workspace
!
   real(r8) zsi(plat)      ! sine of latitudes
   real(r8) zw(plat)       ! Gaussian weights
   real(r8) zra2           ! ra squared
   real(r8) zalp(2*pspt)   ! Legendre function array
   real(r8) zdalp(2*pspt)  ! Derivative array
   real(r8) zslat          ! sin of lat  and cosine of colatitude

   integer i           ! longitude index
   integer j           ! Latitude index
   integer k           ! Level index
   integer kk          ! Level index
   integer kkk         ! Level index
   integer m,lm,mr,lmr ! Indices for legendre array
   integer n           ! Index for legendre array
   integer nkk         ! Print control variables
   integer ik1         ! Print index temporary variable
   integer ik2         ! Print index temporary variable
   integer itmp        ! Dimension of polynomial arrays temporary.
   integer iter        ! Iteration index
   real(r8)    zdt         ! Time step for settau

   logical lprint      ! Debug print flag
   integer irow        ! Latitude pair index
   integer lat         ! Latitude index

   real(r8) xlat           ! Latitude (radians)
   real(r8) pi             ! Mathematical pi (3.14...)
   real(r8) dtime          ! timestep size [seconds]
!
!-----------------------------------------------------------------------
!
   lprint = masterproc .and. .FALSE.

   dtime = get_step_size()
   zdt = dtime/eul_nsplit
   if (masterproc) write(iulog,*) 'zdt, dtime, eul_nsplit=',zdt,dtime,eul_nsplit

!   call hdinti  (rearth  ,dtime   )
   call hdinti  (rearth  ,zdt   )
!
! Initialize commap.  Set hybrid level dependent arrays
!
   call hycoef_init

   if ( .not. single_column ) then
!
! NMAX dependent arrays
!
      if (pmmax.gt.plon/2) then
         call endrun ('INITCOM:mmax=ptrm+1 .gt. plon/2')
      end if
   endif

   zra2 = ra*ra
   do j=2,pnmax
      sq(j)  = j*(j-1)*zra2
      rsq(j) = 1._r8/sq(j)
   end do
   sq(1)  = 0._r8
   rsq(1) = 0._r8
!
! MMAX dependent arrays
!
   do j=1,pmmax
      xm(j) = j-1
   end do
   if ( .not. single_column ) then
!
! Gaussian latitude dependent arrays
!
      call gauaw(zsi     ,zw      ,plat    )
      do irow=1,plat/2
         slat(irow) = zsi(irow)
         w(irow)              = zw(irow)
         w(plat - irow + 1)   = zw(irow)
         cs(irow)  = 1._r8 - zsi(irow)*zsi(irow)
         xlat = asin(slat(irow))
         clat(irow) = -xlat
         clat(plat - irow + 1) = xlat
      end do
      
      do lat=1,plat
         latdeg(lat) = clat(lat)*45._r8/atan(1._r8)
      end do
   endif
!
! Integration matrices of hydrostatic equation(href) and conversion
! term(a).  href computed as in ccm0 but isothermal bottom ecref
! calculated to conserve energy
!
   do k=1,plev
      do kk=1,plev
         href(kk,k) = 0._r8
         ecref(kk,k) = 0._r8
      end do
   end do
!
! Mean atmosphere energy conversion term is consistent with continiuty
! Eq.  In ecref, 1st index = column; 2nd index = row of matrix.
! Mean atmosphere energy conversion term is energy conserving
!
   do k=1,plev
      ecref(k,k) = 0.5_r8/hypm(k) * hypd(k)
      do kk=1,k-1
         ecref(kk,k) = 1._r8/hypm(k) * hypd(kk)
      end do
   end do
!
! Reference hydrostatic integration matrix consistent with conversion
! term for energy conservation.  In href, 1st index = column; 
! 2nd index = row of matrix.
!
   do k = 1,plev
      do kk = k,plev
         href(kk,k) = ecref(k,kk)*hypd(kk)/hypd(k)
      end do
   end do
!
! Print statements
!
   if (lprint) then
      nkk = plev/13
      if (mod(plev,13).ne.0) nkk = nkk + 1
      write(iulog,*)' '
      write(iulog,*)'INITCOM: Hydrostatic matrix href'
      do kk=1,nkk
         ik1 = 1 + (kk-1)*13
         ik2 = min0( ik1+12, plev )
         write(iulog,9920) (k,k=ik1,ik2)
         do kkk=1,plev
            write(iulog,9910) kkk,(href(kkk,k),k=ik1,ik2)
         end do
      end do
      write(iulog,*)' '
      write(iulog,*)'INITCOM: Thermodynamic matrix ecref'
      do kk=1,nkk
         ik1 = 1 + (kk-1)*13
         ik2 = min0( ik1+12, plev )
         write(iulog,9920) (k,k=ik1,ik2)
         do kkk=1,plev
            write(iulog,9910) kkk,(ecref(kkk,k),k=ik1,ik2)
         end do
      end do
   end if
!
! Multiply href by r
!
   do k=1,plev
      do kk=1,plev
         href(kk,k) = href(kk,k)*rair
      end do
   end do
!
! Compute truncation parameters
!
   if (masterproc) then
      write(iulog,9950) ptrm,ptrn,ptrk
   end if
!
! Compute semi-implicit timestep constants (COMSPE)
!
!   zdt = dtime
!   if (nsrest==0) zdt = 0.5_r8*zdt
!
!   now initizized in dynpgk:
!   call settau(zdt)
!
! Determine whether full or reduced grid
!
   fullgrid = .true.
   do j=1,plat
      if (masterproc) then
         write(iulog,*)'nlon(',j,')=',nlon(j),' wnummax(',j,')=',wnummax(j)
      end if
      if (nlon(j).lt.plon) fullgrid = .false.
   end do

   if ( single_column ) then
      do j=1,plat
         slat(j) = 1.0_r8 * sin(4.0_r8*atan(1.0_r8)*scmlat/180._r8)
         w(j)   = 2.0_r8/plat
         cs(j)  = 10._r8 - slat(j)*slat(j)
!         itmp = 2*pspt - 1
!         call phcs  (zalp    ,zdalp   ,itmp    ,zslat    )
!         call reordp(j       ,itmp    ,zalp    ,zdalp   )
      end do

!
! Latitude array (list of latitudes in radians)
!
      xlat = asin(slat(1))
      clat(1) = xlat
      
      clat(1)=scmlat*atan(1._r8)/45._r8
      latdeg(1) = clat(1)*45._r8/atan(1._r8)
      clon(1,1)   = 4.0_r8*atan(1._r8)*mod((scmlon+360._r8),360._r8)/180._r8
      londeg(1,1) = mod((scmlon+360._r8),360._r8)
!
! SCAM not yet able to handle reduced grid.
!
      if (.not. fullgrid) then
         call endrun ('INITCOM: SCAM not yet configured to handle reduced grid')
      end if
   else
!
! Compute constants related to Legendre transforms
! Compute and reorder ALP and DALP
!
      allocate( alp  (pspt,plat/2) )
      allocate( dalp (pspt,plat/2) )
      do j=1,plat/2
         zslat = slat(j)
         itmp = 2*pspt - 1
         call phcs  (zalp    ,zdalp   ,itmp    ,zslat    )
         call reordp(j       ,itmp    ,zalp    ,zdalp   )
      end do
!
! Copy and save local ALP and DALP
!
      allocate( lalp  (lpspt,plat/2) )
      allocate( ldalp (lpspt,plat/2) )
      do j=1,plat/2
         do lm=1,numm(iam)
            m = locm(lm,iam)
            mr = nstart(m)
            lmr = lnstart(lm)
            do n=1,nlen(m)
               lalp(lmr+n,j) = alp(mr+n,j)
               ldalp(lmr+n,j) = dalp(mr+n,j)
            end do
         end do
      end do
!
! Mirror latitudes south of south pole
!
      lat = 1
      do j=j1-2,1,-1
         nlonex(j) = nlon(lat)
         lat = lat + 1
      end do
      nlonex(j1-1) = nlon(1)     ! south pole
!
! Real latitudes
!
      j = j1
      do lat=1,plat
         nlonex(j) = nlon(lat)
         j = j + 1
      end do
      nlonex(j1+plat) = nlon(plat)  ! north pole
!
! Mirror latitudes north of north pole
!
      lat = plat
      do j=j1+plat+1,platd
         nlonex(j) = nlon(lat)
         lat = lat - 1
      end do
!
! Longitude array
!
      pi = 4.0_r8*atan(1.0_r8)
      do lat=1,plat
         do i=1,nlon(lat)
            londeg(i,lat) = (i-1)*360._r8/nlon(lat)
            clon(i,lat)   = (i-1)*2.0_r8*pi/nlon(lat)
         end do
      end do

      do j=1,plat/2
         nmmax(j) = wnummax(j) + 1
      end do
      do m=1,pmmax
         do irow=1,plat/2
            if (nmmax(irow) .ge. m) then
               beglatpair(m) = irow
               goto 10
            end if
         end do
         call endrun ('INITCOM: Should not ever get here')
10       continue
      end do
!
! Set up trigonometric tables for fft
!
      do j=1,plat
         call set99(trig(1,j),ifax(1,j),nlon(j))
      end do
!
! Set flag indicating dynamics grid is now defined.
! NOTE: this ASSUMES initcom is called after spmdinit.  The setting of nlon done here completes
! the definition of the dynamics grid.
!
   endif

   return

9910 format( 1x,i3,13f9.5)
9920 format(/,      13i9)
9950 format(/,'     Truncation Parameters',/,'     NTRM = ',i4,/, &
      '     NTRN = ',i4,/,'     NTRK = ',i4,/)

end subroutine initcom
